• Неметрическое многомерное шкалирование
  • Как работает неметрическое многомерное шкалирование
  • Оценка качества подгонки ординации
  • Трактовка результатов ординации

Вы сможете

  • Построить диаграмму nMDS.
  • Охарактеризовать качество ординации с помощью величины стресса.
  • Сравнить результаты нескольких ординаций.

Проблема изображения многомерных данных

Облако точек в многомерном пространстве

Когда признаков много, можно представить все объекты как облако точек в многомерном пространстве.

Migration by Don McCullough on Flickr

Для изображения N объектов в идеале нужно N-1 измерений

  • 2 объекта = 1D прямая
  • 3 объекта = 2D плоскость
  • 4 объекта = 3D объем
  • N объектов = (N-1)-мерное пространство

Многомерное пространство сложно изобразить. Есть два пути:

  • Выбрать наиболее информативную проекцию (не все можно хорошо спроецировать).
  • Сохранить отношения между объектами (придется исказить расстояния).

Ежик в тумане

black shadows for a white horses / les negres ombres dels cavalls blancs by Ferran Jordà on Flickr

Неметрическое многомерное шкалирование

Nonmetric Multidimensional Scaling (nMDS)

Метод визуализации отношений между объектами в пространстве с небольшим числом измерений (обычно 2).

Исходные данные — матрица расстояний между объектами в многомерном пространстве.

nMDS ординация старается сохранить отношения между объектами.

Карта пригородов Санкт-Петербурга

Расстояния по автодорогам

Расстояния по автородогам не совсем “евклидовы”: например, из Кронштадта до Санкт-Петербурга нельзя добраться по прямой.

Обычная проекция исходного пространства на плоскость неудачна

Карта







Проекция при помощи анализа главных координат (PCoA)

  • Карта неудачна: Стрельна, Красное Село и Гатчина сливаются вместе, хотя они не рядом.

Можно восстановить карту, сохраняя ранги расстояний между объектами

Карта





Ординация nMDS

  • Что странного в этой карте?

Важные свойства nMDS

  1. Ординация сохраняет ранг расстояний между объектами (похожие располагаются близко, непохожие — далеко; если объект А похож на B, больше чем на C, то и на ординации он, скорее всего, окажется ближе к B, чем к C).

  2. На ординации имеет смысл только взаиморасположение объектов. Облако точек в осях MDS можно вращать, перемещать, зеркально отражать. Суть ординации от этого не изменится.

  3. Значения координат точек в ординации лишены смысла (их вообще можно не приводить на итоговой ординации).

Как работает nMDS

Наблюдаемые различия и ожидаемые расстояния

Взаиморасположение точек на плоскости подобно взаиморасположению точек в многомерном пространстве признаков, но значения не полностью совпадают.

Наблюдаемые различия и ожидаемые расстояния

Мы можем построить график зависимости расстояний между точками на ординации и соответствующих им значений коэффициентов различия в исходном пространстве.


  • По оси X — значения коэффициентов различий \(d_{h,i}\), ранжированные в порядке возрастания
  • По оси Y — значения расстояний на плоскости ординации — \(\hat d_{h,i}\)

У хорошей ординации nMDS ранги расстояний должны совпадать с рангами коэффициентов различия в исходном многомерном пространстве.

Диаграмма Шепарда

Диаграмма Шепарда показывает соответствие расстояний на ординации и коэффициентов различия в исходном многомерном пространстве.


  • По оси X — значения коэффициентов различий \(d_{h,i}\), ранжированные в порядке возрастания
  • По оси Y — значения расстояний на плоскости ординации — \(\hat d_{h,i}\)

В идеале, большим расстояниям в исходном пространстве должны соответствовать большие расстояния на ординации, а меньшим — меньшие (т.е. должны сохраняться ранги расстояний).

Красная линия — предсказанные монотонной регрессией значения расстояний на ординации — монотонно возрастает.

Монотонная регрессия

Монотонная регрессия минимизирует сумму квадратов отклонений значений координат на ординации от значений в исходном многомерном пространстве, сохраняя монотонные отношения между ними.

Алгоритм подбора монотонной регрессии называется Pool-Adjacent-Violator-Algorithm (PAVA)

Стресс — мера качества ординации

Стресс (Stress) показывает, насколько соответствуют друг другу взаиморасположение точек на плоскости ординации и в исходном многомерном пространстве признаков.

Стресс — это корень из суммы квадратов остатков монотонной регрессии, отнесенной к сумме квадратов всех коэффициентов различия в исходной матрице.

\[ Stress(1) = \sqrt{\frac{\sum_{h < i}{(d_{h,i} - \hat d_{h,i})^2}}{\sum_{h < i}{d_{{h,i}}^2}}}\]

  • \(d\) — наблюдаемое значение коэффициента различия между двумя точками
  • \(\hat d\) — значение, предсказанное монотонной регрессией

Хорошо, когда \(Stress(1) < 0.2\).

Алгоритм MDS (для двумерного случая)


Ежик в тумане





  1. Вычисляем матрицу коэффициентов различия между объектами.
  2. Распределеяем объекты на плоскости в случайном порядке (или при помощи PCoA).
  3. Вычисляем стресс.
  4. Сдвигаем точки по плоскости так, чтобы минимизировать стресс.
  5. Повторяем шаги 2-4 несколько раз, чтобы избежать локальных минимумов стресса.
  6. Обычно финальную ординацию поворачивают, чтобы вдоль оси X было максимальное варьирование.

Ограничения nMDS

  • nMDS — метод визуализации данных.
  • Коэффициент различия и трансформация данных должны быть выбраны исходя из того, какие именно свойства объектов должны быть визуализированы.
  • Число измерений, выбранное для решения, может влиять на результат.

Пример

Симбионты мидий

Открываем данные

dat <- read.delim("data/Krapivin_2017_Medulis-symb.tab", skip = 36)
str(dat)
## 'data.frame':    548 obs. of  23 variables:
##  $ Event                 : Factor w/ 3 levels "Gorelyi_Isl",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Date.Time             : Factor w/ 1 level "2013-07": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Ord.No                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Area                  : Factor w/ 3 levels "B. Gorelyi isl.",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Site                  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Zone                  : Factor w/ 3 levels "down","middle",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Comment               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ M..edulis.L.shell..mm.: num  31.2 41.1 46.4 40 34.5 47.6 34.3 49.6 49.8 30 ...
##  $ M..edulis.age..a.     : num  1 2 2 3 4 3 2 4 3 1 ...
##  $ Urastoma.spp.....     : int  26 16 19 11 14 16 10 4 11 0 ...
##  $ Renicola.spp.....     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Himasthla.spp.....    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Metacercaria....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gymnophallus.spp..... : int  0 0 0 0 0 0 0 0 4 0 ...
##  $ Bucephalidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alg....               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Nematoda....          : int  1 1 0 0 0 0 0 0 0 0 ...
##  $ Microsetella....      : int  0 0 0 0 0 2 0 0 0 0 ...
##  $ Copepoda....          : int  0 0 0 0 0 1 0 0 0 0 ...
##  $ Chironomidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Halacaridae....       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Jaera.spp.....        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ostrac....            : int  0 0 0 0 0 0 0 0 0 0 ...

Приводим в порядок названия переменных

colnames(dat) # Было
##  [1] "Event"                  "Date.Time"              "Ord.No"                
##  [4] "Area"                   "Site"                   "Zone"                  
##  [7] "Comment"                "M..edulis.L.shell..mm." "M..edulis.age..a."     
## [10] "Urastoma.spp....."      "Renicola.spp....."      "Himasthla.spp....."    
## [13] "Metacercaria...."       "Gymnophallus.spp....."  "Bucephalidae...."      
## [16] "Alg...."                "Nematoda...."           "Microsetella...."      
## [19] "Copepoda...."           "Chironomidae...."       "Halacaridae...."       
## [22] "Jaera.spp....."         "Ostrac...."
colnames(dat) <- gsub("[.]", replacement = "", colnames(dat))
dat <- rename(dat, L = "MedulisLshellmm", Age = "Medulisagea")

colnames(dat) # Стало
##  [1] "Event"           "DateTime"        "OrdNo"           "Area"           
##  [5] "Site"            "Zone"            "Comment"         "L"              
##  [9] "Age"             "Urastomaspp"     "Renicolaspp"     "Himasthlaspp"   
## [13] "Metacercaria"    "Gymnophallusspp" "Bucephalidae"    "Alg"            
## [17] "Nematoda"        "Microsetella"    "Copepoda"        "Chironomidae"   
## [21] "Halacaridae"     "Jaeraspp"        "Ostrac"

Приводим в порядок данные

# Делаем сайт фактором
dat$Site <- factor(dat$Site, levels = c(2, 3, 1), labels = c("G","L","S"))

# Сливаем редкие виды в одну категорию
f_remove <- c("Nematoda", "Microsetella", "Copepoda", 
              "Chironomidae", "Halacaridae", "Jaeraspp",
              "Ostrac")
dat$Other <- rowSums(dat[, f_remove])

# Суммарная численность симбионтов
f_sp <- c("Urastomaspp", "Renicolaspp", "Himasthlaspp",
          "Metacercaria", "Gymnophallusspp", "Alg",
          "Other")
dat$Total  <- rowSums(dat[, f_sp])

# Данные для анализа
# Только мидии с симбионтами и возрастом от 3 до 8 лет 
dat <- dat[dat$Total != 0 & dat$Age %in% 3:8, ]
spec <- dat[, f_sp]                         # виды-симбионты
env <- dat[, c("Zone", "Site", "L", "Age")] # свойства мидий-хозяев

Ординация сообществ симбионтов в мидиях

Задание 1

Постройте ординацию мидий по обилию разных видов-симбионтов с использованием коэффициента различия Брея-Куртиса. Следите, чтобы алгоритму удалось найти финальное решение. Если необходимо, настройте вызов metaMDS()

Вычислите стресс для получившейся ординации.

Нарисуйте простейший график.

ord_mus <- 

Решение

ord_mus <- metaMDS(spec, dist = "bray")

По-умолчанию metaMDS() пытается подобрать оптимальную трансформацию, подходящую для анализа сообществ. В данном случае — извлекает корень и делает двойную висконсинскую стандартизацию.

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1319177 

Такая модель не сходится.

*** No convergence -- monoMDS stopping criteria:
    19: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

Решение

После отключения автоматического подбора трансформации данных удается найти решение.

ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
## Run 0 stress 0.141 
## Run 1 stress 0.145 
## Run 2 stress 0.141 
## ... Procrustes: rmse 0.00388  max resid 0.0344 
## Run 3 stress 0.145 
## Run 4 stress 0.146 
## Run 5 stress 0.144 
## Run 6 stress 0.154 
## Run 7 stress 0.154 
## Run 8 stress 0.141 
## ... Procrustes: rmse 0.00221  max resid 0.0231 
## Run 9 stress 0.146 
## Run 10 stress 0.141 
## ... Procrustes: rmse 0.00273  max resid 0.0447 
## Run 11 stress 0.141 
## ... Procrustes: rmse 0.0034  max resid 0.0459 
## Run 12 stress 0.141 
## ... Procrustes: rmse 0.00311  max resid 0.034 
## Run 13 stress 0.144 
## Run 14 stress 0.141 
## ... Procrustes: rmse 0.00342  max resid 0.0419 
## Run 15 stress 0.141 
## ... Procrustes: rmse 0.00299  max resid 0.038 
## Run 16 stress 0.141 
## ... Procrustes: rmse 0.00318  max resid 0.0343 
## Run 17 stress 0.144 
## Run 18 stress 0.15 
## Run 19 stress 0.141 
## ... Procrustes: rmse 0.00214  max resid 0.0217 
## Run 20 stress 0.141 
## ... Procrustes: rmse 0.00236  max resid 0.0329 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress ratio > sratmax

Решение

Хорошее решение, но график нужно украсить.

ord_mus$stress
## [1] 0.141
ordiplot(ord_mus)

Украшенный график

# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)

# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])

Украшенный график с центроидами видов

ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
text(ord_mus, display = "spec", cex = 0.9, col = "grey20")

Визуализация ординации в ggplot2

Задание 2

Используя данные, приведенные ниже, постройте график ординации при помощи ggplot2. Покажите точками — мидий, текстом — центроиды для видов-симбионтов, цветом — зону литорали, формой — сайт, размером маркеров — размер мидий.

# Координаты точек (мидий)
ord_mus_pt <- data.frame(env, scores(ord_mus, display = "sites"))
head(ord_mus_pt, 2)
# Координаты центроидов переменных (видов-симбионтов)
ord_mus_sp <- data.frame(scores(ord_mus, display = "species"))
ord_mus_sp$Species <- rownames(ord_mus_sp)
head(ord_mus_sp, 2)

Решение: График с точками

gg_ord_mus <- ggplot() +
  geom_point(data = ord_mus_pt, 
             aes(x = NMDS1, y = NMDS2, colour = Zone, 
                 shape = Site, size = L), alpha = 0.5)
gg_ord_mus

Решение: График с центроидами видов

gg_ord_mus_sp <- gg_ord_mus +
  geom_text(data = ord_mus_sp, 
            aes(x = NMDS1, y = NMDS2, 
                label = Species))
gg_ord_mus_sp

Интерпретация ординации: envfit

Анализ связи с внешними переменными c помощью функции envfit()

Функция envfit() строит регрессионную модель зависимости переменной среды от координат точек на плоскости ординации.

Модель вот такого вида:

\(y_i = \beta_1 \text{NMDS1}_i + \beta_2 \text{NMDS2}_i + \varepsilon_i`\)

Значимость зависимости оценивается при помощи пермутационного теста (это позволяет работать в ситуации, когда неполностью выполнены условия применимости).

Такой вызов функции позволит подобрать сразу несколько моделей и изобразить их на графике ординации:

ef <- envfit(ord_mus, env[, c("Zone", "Site", "L", "Age")])

Интерпретация результатов envfit() для непрерывных переменных

ef$vectors
##      NMDS1  NMDS2   r2 Pr(>r)    
## L   -1.000 -0.030 0.19  0.001 ***
## Age  0.790 -0.613 0.04  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999




  • Колонки MDS1 и MDS2 содержат информацию о направлении вектора непрерывной переменной. Длина вектора пропорциональна корню из коэффициента детерминации для модели.
  • r2 — \(R^2\) из регресcионного анализа.
  • Pr(>r) — пермутационная оценка статистической значимости.

Интерпретация результатов envfit() для дискретных переменных

ef$factors
## Centroids:
##            NMDS1 NMDS2
## Zonedown   -0.65 -0.09
## Zonemiddle -0.14 -0.05
## Zoneup      0.83  0.14
## SiteG       0.19  0.17
## SiteL       0.23 -0.41
## SiteS      -0.42  0.22
## 
## Goodness of fit:
##        r2 Pr(>r)    
## Zone 0.37  0.001 ***
## Site 0.16  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
  • Колонки MDS1 и MDS2 содержат координаты центроидов для дискретных факторов (т.е. положение центра облака точек, относящихся к категории фактора).
  • r2 – \(R^2\) из регресcионного анализа.
  • Pr(>r) — пермутационная оценка статистической значимости.

График с векторами и центроидами, найденными envfit()

ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
plot(ef)

Для ggplot-графика удобно использовать вспомогательный пакет

# install.packages("devtools")
# devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
ord_mus_ef <- fortify(ef)
ord_mus_ef

ggplot2 версия графика envfit()

gg_ord_mus +
  geom_segment(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
               aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm"))) +
  geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
            aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1)) +
  geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Centroid", ],
            aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1))

Ограничения envfit()

  • Осторожно! envfit() основан на предположении о линейном характере связи между внешней переменной и координатами на ординации. Это не всегда так.

Интерпретация ординации: ordisurf

Анализ связи с внешними переменными c помощью ordisurf()

ordisurf() применяется, когда зависимость нелинейна.

В основе работы функции — общая аддитивная модель (General Additive Model, GAM).

GAM позволяют подобрать нелинейные зависимости, даже если заранее неизвестна их форма.

ordisurf() использует пакет mgcv для подбора GAM на основе метода тонких пластин (thin plate spline).

График с поверхностью, найденной ordisurf()

Изолинии на графике ординации показывают предсказанное значение зависимой переменной для разных участков ординации.

par(mfrow = c(1, 2))
os_L <- ordisurf(ord_mus, env$L, method = "REML")
os_Age <- ordisurf(ord_mus, env$Age, method = "REML")
par(mfrow = c(1, 1))

Интерпретация результатов ordisurf()

Длина мидий линейно меняется на ординации (эффективное число степеней свободы для сплайна \(edf=1.96\)). Зависимость довольно хорошо выражена (объясненная девианса 18.7%).

summary(os_L)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   43.410      0.484    89.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df    F p-value    
## s(x1,x2) 1.96      9 9.83  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.183   Deviance explained = 18.7%
## -REML = 1457.6  Scale est. = 92.439    n = 395

Интерпретация результатов ordisurf()

Возраст мидий меняется нелинейно вдоль ординационных осей (эффективное число степеней свободы для сплайна \(edf=4.99\)). Эта зависимость довольно слабая (объясненная девианса 8.97%).

summary(os_Age)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8861     0.0713    68.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df    F    p-value    
## s(x1,x2) 5.06      9 3.76 0.00000051 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0791   Deviance explained = 9.09%
## -REML =  704.9  Scale est. = 2.0097    n = 395

Для ggplot-графика понадобится добыть данные о контурах

fortify.ordisurf <- function(model) {
  # Fortifies an object of class ordisurf to produce 
  # a dataframe of contour coordinates 
  # for subsequent plotting in ggplot.
  xy <- expand.grid(x = model$grid$x, y = model$grid$y)
  xyz <- cbind(xy, c(model$grid$z))
  names(xyz) <- c("x", "y", "z")
  return(na.omit(xyz))
}

ord_mus_os <- fortify.ordisurf(os_Age)
head(ord_mus_os, 4)

ggplot2 версия графика с поверхностью, найденной ordisurf()

ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) +
  stat_contour(aes(colour = ..level..)) +
  labs(x = "NMDS1", y = "NMDS2", colour = "Age")

Ограничения ordisurf()

  • Осторожно! Для использования ordisurf() важно, чтобы в ординации участвовало много объектов и они располагались на плоскости ординации более-менее равномерно. Пустоты и отскоки снижают качество предсказаний.

Финальный график

Take-home messages

  • nMDS — один из методов ординации в пространстве сниженной размерности. При ординации этим методом сохраняются ранги расстояний между объектами.
  • Мера качества ординации называется “стресс”.
  • Значения координат не имеют особенного смысла. Имеет значение только взаиморасположение точек.
  • Результат nMDS зависит от выбора меры различия.
  • Существует несколько способов визуализировать на ординации значения внешних переменных. При выборе между envfit() и ordisurf() обращайте внимание на форму связи.

Литература

  • Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
  • Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
  • Oksanen, J., 2015. Multivariate analysis of ecological communities in R: vegan tutorial. http://www.cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf
  • Zuur, A. F., Ieno, E. N., Smith, G. M. Analysing Ecological Data. Springer 2007

Самостоятельная работа

Задание 3

Во всех примерах необходимо:

  1. Разобраться с данными.
  2. Построить ординацию объектов (описаний, проб и т.п.).
  3. Визуализировать связь между полученной ординацией и параметрами среды.
  4. Сделать выводы о наиболее важных факторах.

Источники данных

  1. Обилие полихет-трубкостроителей на литорали Белого моря (данные В.М.Хайтова).
  2. растительные сообщества во франции La Mafragh (данные из работы de Belair et al. 1987, пакет ade4).
  3. Деревья на острове Barro Colorado (данные из работы Condit et al. (2002), пакет vegan).